Kameron Lightheart
MATH 321 Section 2
11/29/18
from matplotlib import pyplot as plt
from scipy.io import wavfile
import numpy as np
import IPython
from scipy import fftpack as fftpack
from scipy.fftpack import fft, ifft
from scipy.signal import fftconvolve
from scipy.fftpack import fft2, ifft2
import imageio
plt.rcParams["figure.dpi"] = 300 # Fix plot quality.
plt.rcParams["figure.figsize"] = (12,3) # Change plot size / aspect (you may adjust this).
class SoundWave(object):
"""A class for working with digital audio signals."""
# Problem 1.1
def __init__(self, rate, samples):
"""Set the SoundWave class attributes.
Parameters:
rate (int): The sample rate of the sound.
samples ((n,) ndarray): NumPy array of samples.
"""
self.rate = rate
self.samples = samples
n = len(self.samples)
indicies = np.arange(n)
self.frequencies = indicies * self.rate / n
# Problems 1.1 and 1.7
def plot(self, myBool=False):
"""Plot the graph of the sound wave (time versus amplitude)."""
n = len(self.samples)
if myBool:
Samples = fftpack.fft(self.samples)
fig, ax1 = plt.subplots(1,2)
# Plot the graph of the sound wave
domain = np.linspace(0, (n/self.rate), n)
ax1[0].plot(domain, self.samples)
ax1[0].set_title("Soundwave Plot")
ax1[0].set_ylim(-32768.0,32767.0)
ax1[0].set_xlabel("Seconds")
ax1[0].set_ylabel("Amplitude")
indicies = np.arange(n)
frequencies = indicies * self.rate / n
self.frequencies = frequencies
magnitude = np.abs(Samples)
self.magnitude = magnitude
ax1[1].plot(frequencies[:n//2], magnitude[:n//2])
ax1[1].set_title("Frequencies vs. Magnitude")
ax1[1].set_xlabel("Frequencies (Hz)")
ax1[1].set_ylabel("Magnitude")
plt.show()
else:
# Plot the graph of the sound wave
domain = np.linspace(0, (n/self.rate), n)
plt.plot(domain, self.samples)
plt.suptitle("Soundwave Plot")
plt.ylim(-32768.0,32767.0)
plt.xlabel("Seconds")
plt.ylabel("Amplitude")
plt.show()
# Problem 1.2
def export(self, filename, force=False):
"""Generate a wav file from the sample rate and samples.
If the array of samples is not of type np.int16, scale it before exporting.
Parameters:
filename (str): The name of the wav file to export the sound to.
"""
self.filename = filename
# Scale samples if not right dtype or force = false
# Then write .wav file
if (force or (self.samples.dtype != np.int16)):
scaled_samples = np.int16((self.samples * 32767.0) / np.max(np.abs(self.samples)))
wavfile.write(filename, self.rate, scaled_samples.real)
else:
wavfile.write(filename, self.rate, self.samples)
# Problem 1.4
def __add__(self, other):
"""Combine the samples from two SoundWave objects.
Parameters:
other (SoundWave): An object containing the samples to add
to the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the combined samples.
Raises:
ValueError: if the two sample arrays are not the same length.
"""
if len(self.samples) != len(other.samples):
raise ValueError("Samples are not the same length!!!")
# Add frequencies together to make chord
new_samples = self.samples + other.samples
return SoundWave(self.rate, new_samples)
# Problem 1.4
def __rshift__(self, other):
"""Concatentate the samples from two SoundWave objects.
Parameters:
other (SoundWave): An object containing the samples to concatenate
to the samples contained in this object.
Raises:
ValueError: if the two sample rates are not equal.
"""
if (self.rate != other.rate):
raise ValueError("Rates don't match!!!")
# Add other note to end of self note frequencies
return SoundWave(self.rate, np.hstack((self.samples,other.samples)))
# Problem 2.1
def __mul__(self, other):
"""Convolve the samples from two SoundWave objects using circular convolution.
Parameters:
other (SoundWave): An object containing the samples to convolve
with the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the convolved samples.
Raises:
ValueError: if the two sample rates are not equal.
"""
A = self.samples
B = other.samples
if (self.rate != other.rate):
raise ValueError("Rates don't match!!!")
# Force A and B to be the same length by adding zeros to end
if (len(A) != len(B)):
if (len(A) > len(B)):
diff = len(A) - len(B)
extra = np.zeros(diff)
B = np.hstack((extra, B))
else:
diff = len(B) - len(A)
extra = np.zeros(diff)
A = np.hstack((A, extra))
# Compute DFT of A and B
fftA = fft(A)
fftB = fft(B)
# Compute hadamard of the two
hadamard = fftA * fftA
# Return inverse DFT of this hadamard
samples = ifft(hadamard)
return SoundWave(self.rate,samples)
# Problem 2.2
def __pow__(self, other):
"""Convolve the samples from two SoundWave objects using linear convolution.
Parameters:
other (SoundWave): An object containing the samples to convolve
with the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the convolved samples.
Raises:
ValueError: if the two sample rates are not equal.
"""
A = self.samples
B = other.samples
if (self.rate != other.rate):
raise ValueError("Rates don't match!!!")
n = len(A)
m = len(B)
size = np.ceil(np.log2(n + m - 1))
diffA = np.int(2**size - len(A))
extraA = np.zeros(diffA)
A = np.hstack((A,extraA))
diffB = np.int(2**size - len(B))
extraB = np.zeros(diffB)
B = np.hstack((B,extraB))
# Compute DFT of A and B
fftA = fft(A)
fftB = fft(B)
# Compute hadamard of the two
hadamard = fftA * fftA
# Return inverse DFT of this hadamard
samples = ifft(hadamard)
return SoundWave(self.rate, samples[:n+m-1])
# Problem 2.4
def clean(self, low_freq, high_freq):
"""Remove a range of frequencies from the samples using the DFT.
Parameters:
low_freq (float): Lower bound of the frequency range to zero out.
high_freq (float): Higher boound of the frequency range to zero out.
"""
# Compute DFT
DFT = fft(self.samples)
n = len(DFT)
# Get low and high indicies to erase
low_freq = np.int(low_freq * n/self.rate)
high_freq = np.int(high_freq * n/self.rate)
# Set values between low and high indicies to 0
DFT[low_freq : high_freq] = 0
# DFT is symmetric so set other ends low and high indicies to 0
DFT[n - high_freq : n - low_freq] = 0
# Inverse DFT back to original samples
IDFT = ifft(DFT)
self.samples = IDFT.real
def play(self, filename=None):
if filename is None:
return IPython.display.Audio(rate=self.rate, data=self.samples)
else:
return IPython.display.Audio(filename)
SoundWave.__init__().SoundWave.plot().scipy.io.wavefile.read() and the SoundWave class to plot tada.wav.# Read in tada.wav file
rate, samples = wavfile.read("tada.wav")
# Create SoundWave object and plot
SW = SoundWave(rate, samples)
SW.plot()
SoundWave.export().export() method to create two new files containing the same sound as tada.wav: one without scaling, and one with scaling (use force=True).IPython.display.Audio() to embed the original and new versions of tada.wav in the notebook.# Play tada normal
IPython.display.Audio("tada.wav")
# Play tada new file, no force
SW.export("newfile1.wav")
SW.play("newfile1.wav")
# PLay tada new file, force
SW.export("newfile3.wav", force=True)
SW.play("newfile3.wav")
generate_note().generate_note() to create an A tone that lasts for two seconds. Embed it in the notebook.def generate_note(frequency, duration):
"""Generate an instance of the SoundWave class corresponding to
the desired soundwave. Uses sample rate of 44100 Hz.
Parameters:
frequency (float): The frequency of the desired sound.
duration (float): The length of the desired sound in seconds.
Returns:
sound (SoundWave): An instance of the SoundWave class.
"""
rate = 44100
num_samples = duration * rate
# Make domain
domain = np.linspace(0, num_samples / rate, num_samples)
# Create samples with sin function
samples = np.sin(2 * np.pi * domain * frequency)
return SoundWave(rate, samples)
#Set frequency variables
A = 440
B = 493.88
C = 523.25
D = 587.33
E = 659.25
F = 698.46
G = 783.99
A_up = 880
# Make 2 second A note
A_note = generate_note(A, 2)
# PLay note
IPython.display.Audio(rate=A_note.rate, data=A_note.samples)
SoundWave.__add__().SoundWave.__rshift__().A_note = generate_note(A, 3)
B_note = generate_note(B, 3)
C_note = generate_note(C, 3)
E_note = generate_note(E, 3)
# Make and play chord
A_C_E = A_note + C_note + E_note
IPython.display.Audio(rate=A_C_E.rate, data=A_C_E.samples)
A_note = generate_note(A, 1)
C_note = generate_note(C, 1)
E_note = generate_note(E, 1)
# PLay A then C then E one second each
A_C_E = A_note >> C_note >> E_note
IPython.display.Audio(rate=A_C_E.rate, data=A_C_E.samples)
simple_dft() with the formula for $c_k$ given below.np.allclose() to check that simple_dft() and scipy.fftpack.fft() give the same result (after scaling).def simple_dft(samples):
"""Compute the DFT of an array of samples.
Parameters:
samples ((n,) ndarray): an array of samples.
Returns:
((n,) ndarray): The DFT of the given array.
"""
n = len(samples)
l_list = np.arange(n).reshape(n,1)
# Build DFT matrix
W = np.exp((-2j * np.pi/n) * (l_list @ l_list.T))
return W @ samples / n
# Make random array into sound
array = np.random.random((128,))
SWarray = SoundWave(rate, array)
# Calculate DFt, compare to Scipy dft
myDFT = simple_dft(SWarray.samples)
correctDFT = fftpack.fft(SWarray.samples)
correctDFT /= len(SWarray.samples)
# Check if close enough to correct
np.allclose(myDFT, correctDFT)
simple_fft().simple_dft(), simple_fft(), and scipy.fftpack.fft().
Print the runtimes of each computation.np.allclose() to check that simple_fft() and scipy.fftpack.fft() give the same result (after scaling).def simple_fft(samples, threshold=1):
"""Compute the DFT using the FFT algorithm.
Parameters:
samples ((n,) ndarray): an array of samples.
threshold (int): when a subarray of samples has fewer
elements than this integer, use simple_dft() to
compute the DFT of that subarray.
Returns:
((n,) ndarray): The DFT of the given array.
"""
n = len(samples) # assumed to be a power of 2
if n <= threshold: #this cutoff to be optimized, also a power of 2
return simple_dft(samples)
else:
f_even = simple_fft(samples[::2]) # FFT of even indexed entries of f
f_odd = simple_fft(samples[1::2]) #FFT of odd indexed entries of f
w = np.exp((-2j * np.pi/n) * np.arange(n))
first_sum = f_even + w[:n//2] * f_odd
second_sum = f_even + w[n//2:] * f_odd
return 0.5 * np.concatenate([first_sum, second_sum])
# Get random 8192x1 array, create sound object
array = np.random.random((8192,))
SWarray = SoundWave(rate, array)
# Time all three functions
%time myDFT = simple_dft(SWarray.samples)
%time myFFT = simple_fft(SWarray.samples)
%time correctDFT = fftpack.fft(SWarray.samples)
correctDFT /= len(SWarray.samples)
# Check if FFT is correct
np.allclose(myFFT, correctDFT)
SoundWave.plot() so that it accepts a boolean. When the boolean is True, take the DFT of the stored samples and plot (in a new subplot) the frequencies present on the $x$-axis and the magnituds of those frequences on the $y$-axis. Only the display the first half of the plot, and adjust the $x$-axis so that it correctly shows the frequencies in Hertz.A_note = generate_note(A,1)
rate = 44100
SW = SoundWave(rate, A_note.samples)
# Plot A and its frequencies vs magnitude plot
SW.plot(True)
SW = SoundWave(rate, A_C_E.samples)
# Plot A minor chord and its frequencies vs magnitude plot
SW.plot(True)
rate, samples = wavfile.read("tada.wav")
SW = SoundWave(rate, samples)
# Plot tada.wav file with its frequencies vs magnitude plot
SW.plot(True)
# Play mystery chord
rate, samples = wavfile.read("mystery_chord.wav")
mystery = SoundWave(rate, samples)
IPython.display.Audio(rate=mystery.rate, data=mystery.samples)
Use the DFT to determine the individual notes that are present in mystery_chord.wav.
# Read mystery file
rate, samples = wavfile.read("mystery_chord.wav")
mystery = SoundWave(rate, samples)
# Plot mystery chord
mystery.plot(True)
# Get mystery chord's magnitudes and indicies
magnitude = mystery.magnitude
frequency = mystery.frequencies
n = len(mystery.samples)
# Sort magnitudes
sort = np.argsort(magnitude[:n//2])[::-1]
# Get top notes
note1 = frequency[sort[0]]
note2 = frequency[sort[1]]
note3 = frequency[sort[2]]
note4 = frequency[sort[3]]
print(note1.real)
print(note2.real)
print(note3.real)
print(note4.real)
#note1 = generate_note(note1.real, 1)
#IPython.display.Audio(note1.rate, note1.samples)
A_note = generate_note(A, 4)
G_note = generate_note(G, 4)
C_note = generate_note(C, 4)
D_note = generate_note(D, 4)
note = A_note + G_note + C_note + D_note
IPython.display.Audio(rate=note.rate, data=note.samples)
The notes are...A C D G !!!
SoundWave.__mul__() for circular convolution.tada.wav.tada.wav and the white noise. Embed the result in the notebook.rate, samples = wavfile.read("tada.wav")
tada = SoundWave(rate, samples)
rate = 22050
white_noise = np.random.randint(-32767, 32767, rate*2, dtype=np.int16)
white_noise = SoundWave(rate, white_noise)
convolved = tada * white_noise
convolved_appended = convolved >> convolved
IPython.display.Audio(rate=convolved_appended.rate, data=convolved_appended.samples)
SoundWave.__pow__() for linear convolution.CGC.wav and GCG.wav using SoundWave.__pow__() and scipy.signal.fftconvolve().SoundWave.__pow__() and scipy.signal.fftconvolve() sound the same.rateCGC, sampsCGC = wavfile.read("CGC.wav")
rateGCG, sampsGCG = wavfile.read("GCG.wav")
CGC = SoundWave(rateCGC, sampsCGC)
GCG = SoundWave(rateGCG, sampsGCG)
IPython.display.Audio(rate=CGC.rate, data=CGC.samples)
IPython.display.Audio(rate=GCG.rate, data=GCG.samples)
# Time both my convolve and scipy's, play convolved using __pow__ method
%time CGCconv = CGC**GCG
%time CGCconvSci = fftconvolve(CGC,GCG)
IPython.display.Audio(rate=CGCconv.rate, data=CGCconv.samples)
Use SoundWave.__pow__() or scipy.signal.fftconvolve() to compute the linear convolution of chopin.wav and balloon.wav.
Embed the two original sounds and their convolution in the notebook.
rate, samples = wavfile.read("chopin.wav")
chopin = SoundWave(rate, samples)
rate, samples = wavfile.read("balloon.wav")
balloon = SoundWave(rate,samples)
conv = chopin**balloon
IPython.display.Audio(rate=rate, data=conv.samples)
SoundWave.clean().noisy1.wav by filtering out frequencies from $1250$-$2600$ Hz. Embed the original and the cleaned versions in the notebook.noisy2.wav. Embed the original and the cleaned versions in the notebook.rate, samples = wavfile.read("noisy1.wav")
noisy1 = SoundWave(rate, samples)
noisy1.clean(1250,2600)
noisy1.plot(True)
IPython.display.Audio(rate=noisyClean.rate, data=noisy1.samples)
rate, samples = wavfile.read("noisy2.wav")
noisy2 = SoundWave(rate, samples)
# Plot
noisy2.plot(True)
# Clean noise
noisy2.clean(1200, 15000)
# Plot cleaned graph
noisy2.plot(True)
# Play audio
IPython.display.Audio(rate=noisyClean.rate, data=noisy2.samples)
vuvuzela.wav by filtering bad frequencies out of the left and right channels individually.rate, samples = wavfile.read("vuvuzela.wav")
vuv1 = SoundWave(rate, samples[:,0])
vuv2 = SoundWave(rate, samples[:,1])
# Clean noise
vuv1.clean(200, 500)
vuv2.clean(200, 500)
# Combine the two clean samples
vuv = SoundWave(vuv1.rate, np.vstack((vuv1.samples,vuv2.samples)))
# Play audio
IPython.display.Audio(rate=vuv.rate, data=vuv.samples)
license_plate.png so that the year printed on the sticker in the bottom right corner of the plate is legible.# read in image
im = imageio.imread("license_plate.png")
plt.imshow(im, cmap="gray")
plt.show()
# Compute DFT
im_dft = fft2(im)
# Fix brightness
im_dft[32:38,98:110] = 99
im_dft[178:184,328:336] = 99
im_dft[110:120,125:140] = 99
im_dft[65:75,195:210] = 99
im_dft[140:155,225:240] = 99
im_dft[100:110,300:310] = 99
# Inverse back to original, now fixed
im_idft = np.real(ifft2(im_dft))
plt.imshow(im_idft, cmap="gray")
plt.show()
The year on the sticker is...13